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We study the thermal conductivity of the one-dimensional Fermi-Hubbard model at finite tem¬ 
perature using a density matrix renormalization group approach. The integrability of this model 
gives rise to ballistic thermal transport. We calculate the temperature dependence of the thermal 
Drude weight at half filling for various interactions strength. The finite-frequency contributions 
originating from the fact that the energy current is not a conserved quantity are investigated as 
well. We report evidence that breaking the integrability through a nearest-neighbor interaction 
leads to vanishing Drude weights and diffusive energy transport. Moreover, we demonstrate that 
energy spreads ballistically in local quenches with initially inhomogeneous energy density profiles in 
the integrable case. We discuss the relevance of our results for thermalization in ultra-cold quantum 
gas experiments and for transport measurements with quasi-one dimensional materials. 


Introduction. Improving our understanding of trans¬ 
port in one-dimensional (ID) strongly correlated sys¬ 
tems (SCS) is an active field in condensed matter theory. 
While in ID, the existence of powerful numerical m 
and analytical methods makes it possible to obtain 
quantitative results (see, e.g., [HISHH])) transport coeffi¬ 
cients are very hard to come by exactly and are challeng¬ 
ing quantities to determine numerically. Early studies 
suggested that integrable systems such as the ID spin-1 /2 
Heisenberg or Fermi-Hubbard model (FHM) may pos¬ 
sess ballistic transport properties at finite temperatures 
[S]. In linear response theory, ballistic dynamics mani¬ 
fests itself through non-zero Drude weights. The so far 
best understood model is the spin-1/2 XXZ chain, for 
which the Drude weight for thermal transport has been 
calculated exactly uniiii], while substantial progress has 
recently been made regarding the spin conductivity m- 

m]. 

The theory of transport in the ID FHM is much less 
advanced and has focussed on spin and charge trans¬ 
port [52H25]. The thermal conductivity can, by using 
the Kubo formula [MIISO], be written as 


ReK(a;) = 27rDth(T)(5(w) -|- Kreg(a;) (1) 


with the thermal Drude weight Dth(T) and a regular 
finite-frequency contribution Kreg(a;). Since in SCS, the 
Wiedemann-Franz law is not necessarily valid, indepen¬ 
dent calculations of charge and thermal transport are re¬ 
quired. 

The formal argument to prove a nonzero Drude weight 
relies on the Mazur inequality 01311131] 


^th > 22^2 1 X 

i 


m 


( 2 ) 


where Ith is the energy-current operator and the Qi are 


local or quasi-local conserved quantities [I11I31I33]- In 
the presence of interactions, nontrivial Qi leading to fi¬ 
nite Drude weights typically exist in integrable models 
0. For instance, for the spin-1/2 XXZ chain, the en¬ 
ergy current /th = Qs itself is conserved (implying that 
Kreg(w) = 0), while for the FHM, /th only has a partial 
overlap with Q 3 (both operators have a similar struc¬ 
ture M) such that while Dth > 0, also K^egi^) ^ 0 
0 . As a consequence, the half-filled FHM realizes an 
unusual behavior: ballistic thermal transport 0 , yet dif¬ 
fusive charge conduction HZlUHl at temperatures T > 0. 

Here, we address the outstanding problem of quan¬ 
titatively calculating ReK(a;) at T > 0 for the FHM 
at half filling by using a finite-temperature density ma¬ 
trix renormalization group (DMRG) method 0 l35l - 
139], previously applied to both charge transport in the 
FHM [57| and transport in quasi-lD spin-1/2 systems 
[H SOI EH- We obtain the energy-current au¬ 

tocorrelation functions Cth{t) = Re {Ith{t)Ith)/L from 
time-dependent simulations. Since C'th(t) saturates fast 
at a time-independent non-zero value, we are able to 
extract (i) the thermal Drude weight and (ii) the reg¬ 
ular part from a Fourier transformation in combination 
with linear prediction [42) . Moreover, we consider the 
extended FHM as an example for a non-integrable model 
and provide evidence that ballistic contributions are ab¬ 
sent, with a diffusive form of the low-frequency Kreg(w). 

The FHM has been realized with ultra-cold quantum 
gases I33H3H- In ultra-cold quantum gases, relaxation 
processes play an important role for reaching thermal 
equilibrium during the state preparation [531 133] and 
thermometry is an open experimental problem |55j . Fur¬ 
thermore, understanding thermalization dynamics and 
non-equilibrium transport as such have been the goal of 
several optical-lattice experiments with Hubbard systems 
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FIG. 1. (Color online) DMRG results for Cth{t) for various 
T > 0 at (a) U/to = 2 and (b) U/tg = 8. Thick lines: V = 0, 
thin lines: V = U/3. The curves for U/to = 8, T/to = 0.66 
are multiplied by a factor of 10. 


[55H55] . We demonstrate that real-space perturbations 
in the energy density spread ballistically in the ID FHM 
at T > 0 while charge diffuses naiiHi, providing a route 
to experimentally observing the qualitative difference be¬ 
tween charge and energy dynamics in this model. 

Definitions. The Hamiltonian of the extended FHM is 
given hy H = hi with local terms 

hi = -toY^ (cLg-Uct + +V{ni- l)ini+i - 1) (3) 

(T 

+ ^ ({nit - \)inii -\) + > 


where cia annihilates a fermion with spin a on site I, and 
nia = U and V denote the onsite and the nearest- 

neighbor Coulomb repulsion, respectively. We use open 
boundary conditions. All results in the main text are for 
half filling n = N/L = 1, where N is the total number of 
fermions. For convenience, we implemented the FHM as 
a two-leg spin-1/2 ladder [5i] . 

We derive the energy current from the continuity equa¬ 
tion [5], leading to Ah = (for the full 

expression, see [31]). At n = 1, particle-hole symmetry 
leads to a vanishing thermopower 131EO] and thus, the 
thermal conductivity stems solely from energy-current 
correlations. 

Numerical method. The thermal Drude weight is re¬ 
lated to the long-time asymptote of the current correla¬ 
tion functions via 


Dth = lim lim 

i—>-oo L—>-oo 


ah(t) 
2T2 ’ 


(4) 


and the regular part of the conductivity defined in Eq. 0 
can be obtained from (C'th(0 = C'th(^) — 2T^_Dth) 


^reg(^) 


1 _ p T 

-——Re/ lim C'th(t). 

Ujj. jq L—yoo 


(5) 


FIG. 2. (Golor online) Thermal Drude weight of the integrable 
model versus temperature for U/to = 0,1, 2,4, 8. Solid line: 
exact result Eq. (j^ for [/ = 0, in excellent agreement with 
DMRG data. Dashed lines: High-T behavior Dth = 
with computed using DMRG. Inset: 27rDth//o at T = oo 
(squares), compared to the lower Mazur bound (solid line). 


Note that the derivation of the Kubo formula for k is 
more subtle than for the charge conductivity (see, e.g., 
Refs. [301 10TH05] and references therein and [31] for a dis¬ 
cussion), while there is also ongoing research on thermal 
and energy transport in open quantum system (see, e.g., 

[HOHHO] ). 

Our finite-T DMRG method, implemented via matrix- 
product states [MZOI, is based on the purification trick 
m (see [33111 [75HZ3 for related work). Thus, we sim¬ 
ulate pure states that live in a Hilbert space spanned by 
the physical and auxiliary (ancilla) degrees of freedom. 
Mixed states are obtained by tracing over the ancillas. 
In order to access time scales as large as possible, we 
employ a finite-temperature disentangler |35j , using that 
purification is not unique to slow down the entanglement 
growth. Moreover, we exploit ‘time-translation invari¬ 
ance’ [37], rewrite (/th(t)7th(0)) = (/th(t/2)/th(-t/2)), 
and carry out two independent calculations for /th(t/2) 
as well as for /th(—1/2). Since the energy current is a 
six-point function, a tMDRG simulation for the energy- 
current autocorrelation is much more demanding than 
it is in the charge case. Our calculations are performed 
with L = 100 sites (see Fig. SI in [34] for an analysis 
of the L-dependence). The ‘finite-time’ error of Kreg(i^) 
can be assessed following m, resulting in the error bars 
shown in the figures. 

Real-time decay of Cth. Typical DMRG results for 
C'th(t) are shown in Figs. 0a) and (b) for U = 2to 
and U = 8to, respectively, and temperatures T = 
oo, 2to, 2to/3. We are able to reach times tto ^ 5. For 
V = Q (thick lines), i.e., in the integrable case, C'th(t) 
rapidly saturates at a constant non-zero value, reflect¬ 
ing the ballistic nature of energy transport in this model. 
The transients are surprisingly short compared to spin 
transport in the spin-1/2 XXZ chain [2] and exhibit os- 
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cillations with a fairly small amplitude. To illustrate the 
behavior in the non-integrable extended FHM, we present 
data for V = U/3 (thin lines), for which our system 
Eq. @ is still in the Mott-insulating phase [75]. The 
DMRG results unveil a much stronger decay of C'th(t) 
compared to the integrable case, yet also longer tran¬ 
sient dynamics before the asymptotic regime is reached 
[see, e.g., the data for T = oo shown in Fig. [^a)]. For 
U/to = 8 and T = oo, the real-time decay of Cth is 
consistent with a vanishing D^^(T), as expected for this 
non-integrable model miiMsi]. 

Thermal Drude weight for V = 0. The fast saturation 
of C'th(t) at a constant and non-zero value allows us to 
extract the temperature dependence of Dt\i{T), displayed 
in Fig. 1^ for U/Iq = 0,1, 2,4,8 (note the log-log scale). 
For U = 0, we compare our data to the exact result [50] 

Ah(T) = y [ekVkf{ek)fe^'‘/'^dk, (6) 

where tk = —2tocos(fc), Vk = —dek/dk, and /(e) = 
1/(1 -I- The agreement is excellent. In general, 

Dtii{T) has a maximum at a t7-dependent temperature 
that shifts to larger temperature as U increases. In the 
high-temperature regime T > to, Dth = (dashed 

lines in the figure), where the prefactor has been 
extracted from the numerical data at /3 = 0. In the sup¬ 
plemental material |34j . we compare the thermal Drude 
weight of the FHM to the one of the Heisenberg chain 
m- The latter describes the low-temperature contribu¬ 
tion of spin excitations to the full Dth of the FHM for 
t/ to at n = 1, while we obtained results for the spin- 
incoherent regime T 3> where charge excitations 

dominate. 

We next study how much of the full spectral weight 
of Re«;(a;) is in the Drude peak by plotting 27rDth/7o 
versus t7/to at /? = 0 in the inset of Fig. where 
Iq = f dojReK(u>). The Drude weight contains the full 
weight Jo only at C/ = 0 and for U/to ^ oo. In the 
former case, this results from the exact conservation of 
the thermal current in the non-interacting case, while 
in the latter case, it is a consequence of a full sup¬ 
pression of any scattering between subspaces with dif¬ 
ferent numbers of doublons as U/to diverges. For a finite 
0 < U/to < oo, 27rDth/Jo < 1 and it takes a minimum 
with 27rDth/Jo ~ 0.92 close to U = 2to, implying that 
at n = 1, the dominant contribution to ReK(a;) always 
comes from the Drude weight. Ref. [9] provides a nonzero 
lower bound for Dth at T = oo by considering only Qo 
(a close relative of Jth EH) in Eq. ([|). By comparison 
to this lower bound (solid line in the inset of Fig.[^, we 
conclude that the position of this minimum can be un¬ 
derstood from the competition of the C/-dependent and 
[/-independent contributions to the Drude weight and to 
the total weight Jq oc {l/h)- Moreover, the lower-bound 
from E] is not exhaustive (see also Fig. S2 in [33] showing 
Dth and the lower bound as a function of n). 



FIG. 3. (Color online) Regular part Kreg(w) for U/to = 4 at 
T = 00 . Main panel: non-integrable case {V = U/3). Dashed 
line: Fit of DMRG data to a Lorentzian. Inset: Integrable 
case {V = 0). The data are shown only for uj/to > 1.4, 
whereas for smaller frequencies, uncertainties become too 
large (see the text). 


Non-integrable model and low-frequency dependence of 
Kreg(w)- Upon breaking integrability, our results for 
C'th(t) indicate a vanishing Drude weight, at least at high 
temperatures and for intermediate values of U/to. This 
raises the question of the functional form of Kreg(<j-’) for 
U 7 / 0. Figure 1^ shows /«reg(w) for U/to = 4 at T = 00 
and V = U/3 ^ain panel) and U = 0 (inset). In the 
non-integrable case, Kreg(w) has a broad peak at zero fre¬ 
quency, which is very close to a Lorentzian (the dashed 
line is a fit to the data). This demonstrates that standard 
diffusion is realized in the extended Hubbard model. In 
the integrable case, we often observe maxima in Ki.eg(w) 
at w > 0, which seem to be related to the charge gap. Due 
to the uncertainties involved in extracting the frequency 
dependence, which are due to the finite times reached in 
the simulations and the extraction of the Drude weight, 
we are not able to resolve the low-frequency regime for 
U = 0. Therefore, the question of whether Kreg(w —0) 
is zero or finite in the integrable case, which has been 
intensely studied for spin transport in the spin-1/2 XXZ 
chain [T1I3I|[53HS3], remains an open problem. 

Spreading of local perturbations. The presence of a bal¬ 
listic contribution in the linear response functions trans¬ 
lates into the ballistic spreading of perturbations in the 
local density [H mi [57HSI]- To illustrate this connec¬ 
tion, we study a Hubbard chain at infinite temperature 
and introduce a perturbation in the local charge density 
at t = O'*". This also causes a perturbation in the energy 
density. We measure the time evolution of both densi¬ 
ties Pch,i{t) = {ni{t)) and pth,i(t) = {hi{t)), presented 
in Figs. Qa) and (b). While the energy density shows 
the typical features of a ballistic dynamics [211 HSj j the 
charge density exhibits a much slower spreading and does 
not form fast ballistic jets, nicely illustrating the different 
nature of energy versus charge transport in this model. 
To become more quantitative, we compute the spatial 
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FIG. 4. (Color online) Local qnenches for f//fo = 8 and 
y = 0. (a-c) At time t = 0, a local excitation (|0) + 11 
) + 14) + 1.1114)) is prepared at two sites in the center of 
an otherwise eqnilibrated system at temperature T = oo. 
Real-time evolntion of (a) the normalized local charge den¬ 
sity [{ni{t)} — l]/no, (b) the normalized local energy den¬ 
sity {hi{t)}/Eo, and (c) the relative probability for double 
occupancy = 2) — l/4]/po. The normalization con¬ 

stants read no = (2 * 1.1^ -|- 2)/(l.l^ -|- 3) — l,ifo ~ 0.025to, 
and po — 1.1^/(1.-I- 3) — 1/4. (d) T-quench {Ti/to = 10, 
T 2 = oo): relative energy density (h;(t))/(—0.245to). 

variances associated with the density pth ,ch,l{t) 

1 ^~" 0 

= T7- ^ 0 )^ kh,ch,z(t) - Pth.ch] (V 

•'vth,ch , 

Z=no 

where Iq is the center of the wave packet, tiq cuts off 
boundary effects, denotes the bulk background 

density, and A/th.ch is the excess particle number or en¬ 
ergy induced by the wave packet. As expected we find 
Scrfi^ = crtij(t) — = 0) oc yet a much slower growth 

for the charge cr^j, oc with 1/2 < a < 1 [see Fig. S4]. 
The determination of the exact exponent would require 
longer times and is related to the low-frequency behavior 
of the charge conductivity, yet clearly, charge dynamics is 
not ballistic. Another illustration for the ballistic energy 
spreading can be obtained in T-quenches |21] . in which 
we embed a region with T 2 into a larger system that is 
at Ti < T 2 , which overall has a homogeneous spin and 
charge density. An example is shown in Fig. [^d) and 
as expected, the variance of this wave-packet grows as 
^crth(t)^ oc illustrating that energy spreads ballisti- 
cally in the temperature quench as well. 

The time evolution of the double occupancy (accessi¬ 
ble in optical-lattice experiments) d(t) = (ni^ni^}(t) is 
shown in Fig. |^c) for the quench of Figs. |^a) and (b). 
The profile exhibits fast ballistic jets and the associated 
variance cr^(t) = ^ E/G-^o)^[(dz(0)-^'"®] (D = '£ 1 (^ 1 )) 
increases approximately quadratically at long times and 
is thus sensitive to the fast spreading of the energy den¬ 
sity. For V ^ 0 (see Fig. S5 for an illustration), the 
variances of both energy and double occupancy increase 


much slower than quadratically in time. As a conse¬ 
quence of the ballistic energy transport in a ID FHM 
we expect the absence of thermalization in related quan¬ 
tum gas experiments. From the long-time behavior of 
the respective width Sath,ch{t), one can extract diffusion 
constants HHun] or Drude weights (see, e.g., [55]) via 
Einstein relations (as we verified for our case), provid¬ 
ing an experimental means of measuring transport coef¬ 
ficients. Such a local, real-space and real-time probe for 
thermal transport has recently been used in experiments 
with low-dimensional quantum magnets [93] . Given that 
a coupling to phonons cannot be avoided in quantum 
magnets [94l [95], quantum-gas microscope experiments 
[46H52j could provide a means of studying energy and 
charge transport in the FHM, which is easier to real¬ 
ize with ultracold quantum gases than the spin-1/2 XXZ 
chain in its massive regime, where a similar coexistence 
of diffusive spin transport |20l m] and ballistic energy 
transport exists |91III1|96|. 

Summary and outlook. We computed the thermal con¬ 
ductivity of the ID FHM using a finite-T DMRG method. 
We confirm the ballistic nature of thermal transport in 
the integrable case and we studied the temperature de¬ 
pendence of the Drude weight. The lower bound for 
Dth from i is not exhaustive implying that more lo¬ 
cal (or even quasi-local m (1111331) conserved quan¬ 
tities than just Q 3 play a role. We further demon¬ 
strated that the coexistence of diffusive charge transport 
and ballistic thermal transport is directly reflected in lo¬ 
cal quantum quench dynamics, presumably accessible to 
fermionic quantum gas microscopes [46H52j . For the ex¬ 
tended Hubbard model, we identified regimes in which 
first, the Drude weight clearly vanishes as system size 
increases and second, the low-frequency dependence is 
compatible with diffusive dynamics. 

From the theoretical point of view, an exact calcu¬ 
lation of k{uj) exploiting the integrability of the model 
constitutes an open problem. In view of the existence of 
quasi-ID materials described by the (extended) Hubbard 
model (including some Beechgard salts (97] |98| , organic 
materials [99lH01j . anorganic systems such as Sr 2 Gu 03 
m and carbon nanotubes [103H105] ). a detailed anal¬ 
ysis of energy transport and the calculation of diffusion 
constants is desirable. Finally, the investigation of ther¬ 
moelectric effects in SGS are a timely topic (see, e.g., 
[106H111] ) and should be feasible with our technique, at 
least at high temperatures. 
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FIG. SI. (Color online) Energy-current autocorrelation func¬ 
tion at U/to = 8 and V — 0 calculated for various system 
sizes L. Note that the curves for L = 40 and L = 100 are 
indistinguishable. 


FIG. S2. (Color online) Filling-dependence of the Drude 
weight at U/to = 2 and T = oo. The lower Mazur bound 
from Ref. El is shown as a comparison. 


SI. Energy current operator in fermionic and spin 
language 


Following Ref. O we define the energy current op¬ 
erator using a continuity equation. This yields Ith = 
For V = 0, this results in 


-^th = X! ^0 

l,(T 


(*cj’+i^C;_i^-bh.c.) 

Ur- ■ \r Itl 

2 — + Jlcr) [nicr , 


(SI) 


where jia = ttoC;+icrGcr + local charge current. 

Ith has the same structure as the conserved charge Q 3 , 
except for a different prefactor of 1/2 in front of the U- 
dependent term [S]. 

Within the DMRG numerics, we implement a spin ver¬ 
sion of the Hamiltonian as well as of Ith, which we obtain 
via a Jordan-Wigner transformation: 


Qt = K---ar_i)ar, C4=(naf)(r/...r/_,)rr , 

(S2) 

where and denote standard Pauli matrices, 

and af = {af ± iaf)/2, rf" = {t// ± iT'j')/2. This leads 
to (for notational simplicity we again set P = 0 ) 


H = 


i 


[to(c 


•+ 


U 


i+t(Ji +h.c.)-b-crfr/ 


(S3) 


and 


+ h.c.) 

I 

Oi-i + + {ji-i + ji)^i 1 1 


u r 
4 


(S4) 


where // = —-b h.c., jj = —itoTp/iT; + h.c.. 


-I- 


S2. Transport coefficients 

In general, the heat current Iq of an electronic system 
is defined as [291 HE] 

Iq=Ith-Ul- (S5) 

As a consequence, one needs to consider a two-by-two 
matrix of transport coefficients that relate external 
forces to currents. One possible choice is to work with the 
heat current Iq and the particle(electrical) current I and 
the corresponding forces Fj = — AV (4 -b P), where P 
is the electrostatic potential, 4 is the chemical potential 
and Fq = V(l/T), resulting in 



(we set the electrical charge to one for simplicity). The 
transport coefficients L 4 are computed from the respec¬ 
tive Kubo formulae [29]. 

The thermal conductivity k is measured under the con¬ 
dition of (/) = 0 and, therefore, is given by 

K = L 22 — L2iLti Li2 ■ (S7) 

Alternatively, one can work with the currents Ah and I 
and corresponding transport coefficients Mij. The ther¬ 
mal conductivity is then given by 


K — M22 ~ M 2 lMt^M i 2 . 


(S 8 ) 


Thus, in principle, two terms contribute to k, the second 
arises because of the thermoelectric coupling. As is dis¬ 


cussed in textbooks, Eqs. (S7) and (S 8 ) yield equivalent 
results for k. 

Let us hrst consider the corresponding Drude weights 
Dij associated with the Mij (in the notation of the main 
text. Dll = Ds and D 22 = Dth)- In the case of a 
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FIG. S3. (Color online) Drude weight of the Hubbard model 
at large U/to = 8 in comparison with the result for an 
isotropic Heisenberg spin chain with an antiferromagnetic 
coupling J = 4to/t/. 



FIG. S4. (Color online) Variances ch,d(f) for flis data 
shown in Figs. 4(a-c) of the main text. The dashed lines 
show fits to 5a1{t) ~ {y = th,ch,d) for times tto > 3. The 
fitted exponents read Oth = 1.96, Och = 1.28, and Qd = 1.76. 


half-filled system, e.g., the Mott insulator in ID that 
is considered in the main text, the energy and parti¬ 
cle current have a different parity under particle-hole 
symmetry and therefore, the current-correlation function 
vanishes identically. As a conse¬ 
quence, Mi 2 = M 21 = 0 at half filling at all temper¬ 
atures. This result is known from studies of, e.g., the 
thermopower which consequently also vanishes in a Mott 
insulator at n = 1 in the presence of particle-hole sym¬ 
metry iBOlfTOB]. In other words, since we wrote down the 
Hamiltonian (3) in an explicitly particle-hole symmetric 
form for V = 0, the chemical potential vanishes at half 
filling. For the purpose of our work we conclude that the 
Drude weight Dq associated with the heat current Iq is 
given by 


Dq = Ah . (S9) 

Therefore, for the case studied in the main text, 
namely half filling, the thermal conductivity solely stems 
from the energy-current correlation function Cth(^) = 
(/th(t)/th). 

Note that the justification of the Kubo formula for the 


thermal conductivity |30j is more subtle than for charge 
transport (see [65) for a comprehensive discussion and ad¬ 
ditional references). The main reason is that a gradient in 
temperature is not a mechanical force in the same sense 
as a voltage or gradient in chemical potential is. There¬ 
fore, one cannot add an additional term to the Hamilto¬ 
nian that would contain the temperature gradient, which 
one would commonly treat as the perturbation driving 
the system out of equilibrium. An approach to circum¬ 
vent that problem is to introduce so-called polarization 
operators (see, e.g., [M] for a discussion). Such polariza¬ 
tion operators (essentially of the form P = — can 

then take the role of the external force. This construc¬ 
tion, however, requires open boundary conditions. Other 
attempts make an assumption about local equilibrium or 
rely on the so-called entropy production argument (see 
again [5S] for details). 

Ref. [HS] provides a derivation of the Kubo formula 
that does not rely on the assumption of local equilibrium. 
That paper also attempts a direct comparison between 
time-dependent numerical solutions of certain models in 
the presence of a temperature gradient and mostly con¬ 
firms the validity of the Kubo formula. Some discrep¬ 
ancies were observed in Ref. [65], which, however, were 
later understood to be finite-size effects |113j . 

There is also substantial interest in extracting ther¬ 
mal currents and temperature profiles from simulations 
of open quantum systems (see, e.g., [SSHMl lll4H117j ). 
where typically, master equations are solved. In some 
of these studies (see, e.g., [nsi), however, transport is 
boundary-AvYveT^^ which realizes a different physical sit¬ 
uation from what underlies the derivation of Kuba for¬ 
mulae for charge, spin and energy transport, where the 
bulk experiences a gradient of some force. As a conse¬ 
quence, such open quantum system simulations can yield 
qualitatively different results compared to the Kubo for¬ 
mula, even for charge or spin transport (compare, e.g., 
pn 1116] ). Moreover, it is important to keep in mind that 
conductivities are bulk properties by definition. Thus, 
any contact regions must be excluded in measuring cur¬ 
rents and, in particular, temperature differences. There¬ 
fore, experiments commonly use a four-terminal setup to 
measure the thermal conductivity (see, e.g., [mHHO]). 
The technical challenge for such open quantum system 
simulations is therefore to be able to reach sufficiently 
large systems, as a result of which a direct comparison 
to the Kubo formula is often not possible. 


S3. Comparison of DMRG results for different 
system sizes 


In Fig. SI we present DMRG results for Ah(t) for sev¬ 
eral different system sizes and temperatures at U/t^ = 8 
and V = 0. Since we are using open boundary conditions, 
finite-size effects manifest themselves by a drop of Cth{t) 
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to zero and a sign change at a system-size dependent time 
t* (L). Very importantly, for t < t* (L), the data from dif¬ 
ferent system sizes all coincide, implying that the results 
are for the thermodynamic limit. Moreover, using data 
for open boundary conditions, one needs at least L > 20 
to get to times of the order of 4 ~ I/to. Using L = 100, 
we clearly never reach t *, hence the time-dependent data 
from such a large system is free of finite-size effects in the 
time window reached in the simulations. Note that the 
numerical effort is essentially independent of L but pri¬ 
marily depends on the entanglement growth that limits 
the accessible times. 


S4. Comparison of the Drude weights of the FHM 
and Heisenberg chain 


In principle, two types of excitations should contribute 
to the thermal conductivity. First, for temperatures 
T > AMott, where AMott is the Mott gap, charge ex¬ 
citations (i.e., excitations in the upper Hubbard band) 
become relevant. Second, at very low T <C AMott, charge 
excitations are frozen out and then gapless spin excita¬ 
tions should be the only available ones. Whether these 
two regimes can be resolved depends on the ratio of the 
characteristic energy scale J for magnetic excitations and 
the Mott gap Amom- In the large U/to limit, J = '^tQ/U, 
much smaller than AMott U. Therefore, the very low 
temperature dependence of Dth of the Hubbard model 
should be identical to the one of the spin-1/2 Heisen¬ 
berg model [iniins], where Hth increases linearly at low 
T, takes a maximum at T ;< J/2 and then decreases 
to zero with 1/T^. As a consequence, the full Dth for 
large U/to could have a double maximum structure and 
moreover, at T ^ J, the conductivity should solely stem 
from charge excitations since the spin system will effec¬ 
tively be at T/J = oo (This is often referred to as the 
spin-incoherent regime, see, e.g., [12111122| l. The Drude 
weights of the FHM at U/to = 8 and the spin-1/2 Heisen¬ 
berg chain [lOj are shown in Fig. S3 confirming the qual¬ 


itative picture described above. Note also the difference 
in the maximum values of the Drude weights: the charge 
dominated Dth exceeds the Heisenberg contribution sig¬ 
nificantly. The problem with numerically resolving the 
interesting crossover regime is that it requires a large U, 
which in the units of the FHM renders J a very small 
temperature that at present cannot be reached. 
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